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ABSTRACT 


A computational simulation of the opposed-jet diffusion flame is performed to 
study its structure and extinction limits. The present analysis concentrates on the 
nitrogen-diluted hydrogen-air diffusion flame, which provides the basic information for 
many vehicle designs such as the aerospace plane for which hydrogen is a candidate 
as the fuel. The computer program uses the time-marching technique to solve the 
energy and species equations coupled with the momentum equation solved by the 
collocation method. The procedure is implemented in two stages. In the first stage, a 
one-step forward overal chemical reaction is chosen with the gas phase chemical 
reaction rate determined by comparison with experimental data. In the second stage, 
a complete chemical reaction mechanism is introduced with detailed thermodynamic 
and transport property calculations. Comparison between experimental extinction 
data and theoretical predictions is discussed. The effects of thermal diffusion as well 
as Lewis number and Prandtl number variations on the diffusion flame are also 
presented. 
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I. INTRODUCTION 


The study of stationary flames, as opposed to propagating or explosive flames, 
may be generally divided into two classes, the premixed flame and the diffusion flame. 
In the premixed flame, the fuel and air or oxygen are premixed before entering the 
reaction zone. Normally, the gases of premixed flames start with a high chemical 
potential, which breaks down suddenly in the reaction zone. When the composition 
(the kinds of fuel and oxidizer and their mixture ratio) and the physical conditions (the 
pressure and the temperature) are specified, the final combustion state and the 
characteristics of the flame can be determined uniquely. If mixing occurs rapidly 
compared with combustion reactions or well ahead of the flame zone, burning can be 
considered in terms of homogeneous process, or premixed flame. However, there are 
systems in which mixing controls the burning rate. Most practical systems fall in this 
category, and they are the so-called diffusion flames. 

In the diffusion flame, the fuel and oxidant are initially separated and the 
reactants mix in the same region in which reactions take place. For the diffusion flame, 
unlike the premixed flame, we have no high chemical potential (there may be a small 
potential on the fuel side due to a delay in its breakdown into its elements). Therefore, 
diffusion flames differ from premixed flames in that the combustion occurs at the 
interface between the fuel and the oxidizer, and the burning processes are more 
dependent on the rate of inter-diffusion of the fuel and oxidant than on the rates of the 
chemical reactions involved. For this reason, the aerodynamic nature of the flow is an 
important factor in the behavior of the diffusion flame. The main difficulty with 
diffusion flames is that they do not propagate and, therefore, there are no fundamental 
characteristics like flame velocity which can be readily measured in premixed flames. 
This is also part of the reason that diffusion flames have received less attention than 
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premixed flames in earlier years, despite the fact that diffusion flames have greater 
practical applications and are encountered more frequently. 

The problem of ignition and extinction of flames is of considerable importance to 
re-entry vehicles; the development of flame due to boundary layer heating of these 
vehicles is undesirable and in some cases dangerous. If one could predict the extinction 
speeds for flames in such systems, it would, perhaps, be possible to keep the missile in 
a velocity spectrum beyond the extinction speeds and thus prevent the onset of flames. 
Studies of extinction of laminar counterflow diffusion flames will provide a simple 
model for an accurate understanding of various combustion characteristics closely 
related to the flame's extinction phenomena as well as a detailed understanding of the 
structure of the flame. Laminar counterflow diffusion flames can be classified into two 
large groups : (1) Diffusion flames formed by fuel and oxidant coming from two 
opposed jets (as shown in Fig. 1). and (2) Diffusion flame in the forward stagnation 
region of a porous burner, in which a fuel gas is blown from the porous burner into 
an oncoming oxidizer flow (as illustrated in Fig. 2). 

In the present analysis, only the opposed-jet counterflow diffusion flame 
configuration is selected to analyze the nitrogen-diluted hydrogen-air diffusion flame. 
However, it is easy to apply the present method to the porous burner case with a few 
simple modifications. The boundary layer similarity solution is employed with the use 
of time-marching finite difference technique and one a step second order reaction 
model. The strain induced extinction at high velocity gradient is investigated and 
compared with experimental data. The effects of transport property variations on the 


flame are also examined. 
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Figure 2. Flow Configuration for the Porous Burner. 
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II. LITERATURE SURVEY 

The main purpose of the present theoretical analysis of the diffusion flame is to 
examine how the flame zone thickness and the flame temperature vary with the strain 
induced velocity gradient and/or the reaction rates to determine the critical values of 
the parameters for flame extinction ; it is also the purpose of the present study to 
understand clearly the structure of the diffusion flame. Studies on opposed-jet 
diffusion will help in the evaluation of the chemical kinetic data for fuel-air mixture, 
especially those which are difficult to handle in the premixed flame. 

Some of the pioneering experimental studies on opposed-jet diffusion flame were 
employed by Potter and Butler [1] and Potter, Heimel, and Butler [2| to obtain 
information on overall burning rates of a wide variety of fuels. Spalding [3] made 
analytical studies and obtained an approximate expression for the mass rates at 
extinction. In his model, he assumed infinitely fast reactions to find the location of the 
flame front, the rate at which these reactants are consumed, and the volumetric heat 
release rates. According his theory, the pressure and jet diameter will affect the 
extinction limits of the flame and he recommend that higher Reynolds number should 
be used in the experiment to ensure the diffusion flame is well inside a region one-tenth 
the jet diameter in size, so that the flow can be predicted using his theory. The use of 
variable properties and more complex chemical kinetic models were suggested. Some 
of the conclusions of his theory were later verified by the experimental work of 
Anagnostu and Potter [4]. 

Fendell [5] considered a model for the extinction of the opposed jet diffusion 
flame. He used a direct one-step chemical kinetics of finite rate and very high activation 
energy. The last assumption made it possible to employ asymptotic expansions in his 
analysis. Chung, Fendell and Holt [6] used a one-step reversible model and found the 
conditions which distinguish the simple and multiple transition between frozen and 
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equilibrium state for a given fuel-oxidant system. By employing the asymptotic method, 
they obtained the various ignition and extinction limits governed by the combination 
of Damkohler number and equilibrium constant. 

Jain and Mukunda [7] considered a compressible flow with Prandtl and Schmidt 
number taken as unity. By introducing a conserved property parameter, they reduced 
the range of integration to simplify the numerical analysis. Results for the effects of jet 
temperature, activation energy and oxidant concentration on the extinction condition 
were obtained. Jain and Mukunda [8] later considered a system of hydrogen-carbon 
monoxide fuel mixture. Two one-step reaction schemes were considered and the study 
revealed that a simple mixing rule can be employed to obtain the mass flow rates at 
extinction. The computed maximum volumetric heat release rates were compared with 
experimental data. 

In their work on the prediction of flame location in the stagnation point flow with 
hydrogen injected form porous wall, Liu and Libby [9,10] considered the role of 
variable density through the use of a similarity transformation. Both an infinitely fast 
reaction [9] and a detailed reaction set [10] were used to describe the oxidation of 
hydrogen to water. A single diffusion coefficient was assumed to reduce the numerical 
complexities. In the first case, the boundary layer was considered as an inner zone and 
an outer zone, where oxidant and fuel mass fractions 

were close to zero and were contiguous at the flame sheet. For large rates of injection, 
the governing equations, which were in terms of flame sheet location, were solved by 
the application of asymptotic expansions with different boundary conditions for two 
different zones. The effective heat of formation of the product (H 2 0) was modified to 
account for the creation of secondary products. The results showed good agreement 
with full equilibrium calculation. In the latter case, a quasilinearization was applied 
to solve the stiff boundary value-type differential equations and the authors reported 
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encountering numerical difficulties as the finite-rate chemical reaction rates increased. 
The limiting cases of frozen and equilibrium behavior are discussed. 

Combustion and extinction phenomena in the forward stagnation region of a 
condensed fuel was studied experimentally and theoretically by Tien et al. [1 1]. In their 
numerical analysis, they assumed a second-order forward overall chemical reaction in 
gas phase, with gas-phase activation energy and modified frequency factor, determined 
by comparison with experimental results. The effect of external radiation on the 
extinction limit was also considered. Favorable agreement between experimental 
extinction data and theoretical prediction was obtained for a specified activiation 
energy and modified frequency factor. 

Hahn, Wendt and Tyson [12] used a similar boundary layer approach to get the 
form of the equations in cylindrical coordinates. The conservation equations 
describing a flat laminar opposed-jet, moist CO, diffusion flame was solved by 
numerical integration. Prediction utilizing finite-rate, elementary combustion kinetics 
gave good agreement with a one-dimensional flame experiment. 

Pellett et al. [18,19,20] conducted a series of opposed jet burner experiments to 
study the effect of air-contaminants on hydrogen diffusion flames. They observed a 
new phenomenon, which they named RESTORE. RESTORE is the opposite of 
extinction (or BLOWOFF) in the sense that the flame reestablishes itself at the axis. 
Their results show that the N 2 -contaminated fuel has a more significant effect on the 
BLOWOFF than on the RESTORE condition and, like CO or C0 2 , each has 
moderately small but not negligible thermodynamic and chemical kinetic effects on the 
hydrogen diffusion flame in air. For the water-contaminated hydrogen diffusion flames, 
the results indicate that water will decrease the maximum sustainable H 2 mass flux, just 
prior to extinction, of a N 2 -diluted, H 2 -air diffusion flame. 
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The structure of several methane-air and hydrogen-nitrogen-oxygen counterflow 
diffusion flames were computed by Dixon- Lewis et al. [13]. They also used the 
boundary layer similarity solution and employed complex chemical reactions and 
detailed transport property calculations. Their effort is known to be the first to include 
detailed transport property calculations. They solved the momentum equation by a 
damped Newton method which runs in parallel with a time-dependent approach for the 
coupled energy and species equations. The extinction limits of both the methane-air 
and hydrogen-air flames at high velocity gradient were predicted and the freezing of the 
oxidant reactions on the fuel side of the hydrocarbon flames is discussed in terms of 
their chemical mechanism. Their results of methane-air flame agree well with the 
experimental data. However, the hydrogen-air flame experimental data available at that 
time (Pellett et al. [21]) were overlooked by the authors and, therefore, no comparisons 
were made with experiments. 
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III. MATHEMATICAL MODEL 

The opposed-jet diffusion flame is considered to be an important configuration 
to obtain the kinetic data for highly reactive fuel-oxidant mixtures. Fig. 1 shows the 
geometry considered in the present analysis. We assume the flow is laminar everywhere. 
For nonreacting flow, the solution to this type of flow is given by the solution of two 
inviscid incompressible jets impinging on one another. However, the fuel contained in 
one jet will diffuse into the other jet and the oxidant will do likewise in the opposite 
direction, thus, establishing a narrow zone in which fuel and oxidant coexist. If the 
mixture is ignited, a thin flame surface will be established. The temperature will rise 
sharply at the flame. This narrow temperature peak will affect the density of the gases 
and the flow characteristics will now be different from those of the nonreacting flow. 
This flow must be determined from the solution of the equation of motion coupled with 
the equation of conservation of energy and species. 

Because the effect of the flame on the flow field is restricted to a very narrow 
region, the problem can be considered as being similar to a boundary layer type flow, 
in which the viscous effects are restricted to a small region when compared to the entire 
flow field. This analogy, together with the assumption that the jets are infinitely wide 
compared with the zone of interest, suggests that equations representing the opposed 
jet diffusion flame may be converted into a set of ordinary differential equations with 
the use of a similarity transformation. 

A. GOVERNING EQUATIONS 

A complete description of system shown in Fig. 1 involves the simultaneous 
solution of the equations of conservation of momentum, energy and individual species. 


The general equation of motion is 
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P~jyf = -gradp + div {p grad q) + B + 'F (1) 

where p is the mixture density, D the substantial derivative, q the mass average 
velocity, t the time, p the pressure, p the viscosity, B the body force, and Y stands for 
the viscous terms that are in addition to those expressed by div (p grad q). 

By applying the boundary-layer approach, the corresponding steady state 
momentum equation becomes 

puu x + pvu y = -p x + ( pu y ) y (2) 

where u and v are the velocity components in the x and y directions and the subscripts 
x and y denote differentiation with respect to those independent variables. 


The energy equation may be written for a moving element of fluid as 


P qgrad - div {k grad T) - div QTp q t c t 

— + p div q + O 


where c„ e„ h„ hi are, respectively, the mass fraction, the specific internal energy, 
enthalpy, and enthalpy of formation of species /, k is the thermal conductivity, T is the 
temperature, vv/" is the mass production rate of species i per unit volume, q, is the 
diffusion velocity of species i with respect to q, and <I> is the dissipation function. 

Also, the species conservation equation for species i is given by / 

div{p(q + = w(" (4) 


By combining Equations (3) and (4) with the auxiliary 
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relations 


5>- 


(5) 


Yj^ C ‘ = 0 


(6) 


Z ^ = e ‘ + ^< r 


(7) 


where R, is the gas constant for species i . 


The steady state energy equation reduces to 


pq grad jZ c, (/i, + ^ )j- = div j k grad T - Z pqfi {h t + /i, )j 

+ q grad p + O 


( 8 ) 


With the general boundary layer assumption and expressions for diffusional 
velocity and mixture enthalpy given below 


q, = ~ (A/ c ; ) S rad c i ~ (A 7 )T) Sfad T 


(9) 


h = Z c< (^ + h ‘) 


(10) 


where D , and D, T are, respectively, the molecular and thermal diffusion coefficient of 
species i , we see that Equation (8) becomes 


puh x + pvh y = [KT y ) y + up x + p{u y f + |ZAP { h , + h i) c iy\ y 
+ ( h i + h 'i) T yl T } 


( 11 ) 
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Noting that h, is a function of temperature alone, we can have following 
expression from Equation (10) 

grad h = q {dhJdT ) J grad T + ^{ h i + h i ) Z rad c, (12) 

Let 

r v = X ] C i( dfl il dT ) = X C ' C P ' 1 ( I3 ) 

We can combine Equation (2), (11) and (13) to obtain the desired steady state 
energy equation 


K* + 4 ), + K* + -t\ 

= [(«lc^)(h + u 2 /2)J + {(1/2)Cm - (K/cp](« 2 )3 

rxn . , (14) 

+ {2/ap - (k/c p )ii(^ + ^ ^ j 

+ {^(A^c,/7')(A i + h])T^ 

For the species equation (4), we apply a similar boundary layer analysis and get 
puc ix + pvc iy = { Dfic iy + Dfpc t T y jT + w{" (15) 

For the velocity, v in the>> direction, we introduce the general continuity equation 

C oux\ + (p V x n ) y = 0 (16) 


where n = 0 for two dimensional flow and 1 for axisymmetric flows. 



13 


Equations (2), (14), (15) and (16) thus constitute the system whose solution is 
required. 


B. SIMILARITY TRANSFORMATION 


As is usual in boundary layer problems, we seek solution of restricted form which 
permit reducing the original partial differential equations to ordinary differential form. 
In this problem, we assume that the boundary and initial conditions and the chemical 
behavior are such that velocity and enthalpy profiles remain similar to themselves 
along the stagnation plane, i.e., that all flow and gas parameters are functions of one 
coordinate rj after the coordinate transformation. The procedure is as follows. 


Let 


V = 


(2£)' /2 1 Pe y 


(17) 


and 


-r 


In i 

PefleUeX dx 


(18) 


where subscript e denotes the the edge of boundary layer. 

If we substitute the potential flow solution near the stagnation point (i.e. 
u, — ajt) into the transformation, we get 


V = 


p e a{n + 1) 


m 


dy 


(19) 


14 


and 


x 2(n+,) 

* = PeMe ° 2(« + 1) 


( 20 ) 


The coefficient a in the potential flow solution gives a measure of the 
characteristic time of the flow (Note that it has the dimension, 1/sec) and it depends 
on the flow geometry. The values of a for two dimensional and axisymmetric flow are 
[7] 

a = uj d , for axisymmetric flow 


a — n x uj 4 H , for two dimensional flow 

where d = diameter of the jet and H = width of the jet. 

The transformation from x,y coordinates to £,r] coordinates is carried out by 
means of the relations 


8 _ d 

fy ~ ( 2^) 1/2 dr i 


_d_ 

dx 


2 n 

Pede^e* 




dt 


+ 


d d >1 " 

dn H _ 


( 21 ) 


( 22 ) 


Introducing the stream function 1 1/(£, q) 


+ = (2 


(23) 


the overall continuity equation (16) is automatically satisfied by the usual 


relations 
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dif/ 

pux = IT 


(24) 


n 

pvx 


dji_ 

bx 


(25) 


In addition to Equations (21) — (24) and u/u, = /'(*/) . we choose the following 
dimensionless dependent variables 

* = £ (26, 

r,- 4- (27) 

% 

9 = -L- (28) 

* e 

The assumption that the kinetic energy u 2 l 2 « h 

has been made and this term has been dropped from the energy equation. 

By substituting into Equations (2), (14), (15) and (16) and removing the 
insignificant terms [Appendix A], we have the following similarity equations 

)+//„ + -^(-TT -/,’) = 0 (29) 



+ f St, + 



d_e 
9 drj 



= 0 (30) 
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r 

— L Y- 

p ^ e tin 
1 r 


+ fYn + 
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L t Y- 

L-e, 1 1 


dd 


P r 9 dr\ 


+ 


1 


Wi 


(n+ l)a P c i, 


= 0 


(31) 


pv = - PePeUgX"/ 1 yj 2£ (32) 

where C = pp\p t p„ P, is the Prandtl number {c f pjk), L ti and are the Lewis number 
(D,Pc>/k) and the thermal Lewis number (DfpFJk), respectively. 

These equations, Equations (29) through (31), will be solved by the collocation 
method and time-marching finite difference technique. Equation (32) will be used to 
find normal velocity v. 


C. BOUNDARY CONDITIONS 


1. Momentum Equation . The momentum balance condition will be applied in 
the stagnation region, i.e. 


P- oo U -o 


= Poo U oo 


Therefore 


U - CO) 



Other boundary conditions will be 


/,( oo) = 1 


m = o 


(33) 


(34) 


2. Energy Equation . 
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*(°°) = ^jSi,oo Y i, 


(35) 


3. Species Equations . 


y ; (-oo) = Y,^ , Y t { oo) = r f> 


(36) 


D. THERMODYNAMIC AND TRANSPORT PROPERTIES 

For the present analysis, the thermodynamic and transport property calculations 
are performed by the implementation of the CHEMKIN routines [14] and its extension 
[15]. The CHEMKIN routines provide the thermodynamic properties for the chemical 
species in terms of polynomial fits to the specific heat at constant pressure. Other 
thermodynamic properties are then given in terms of the fits to c p . In the CHEMKIN 
routines, seven coefficients for two temperature ranges are used. These fits are given in 
the following form 


R 


= A] ,• + ^2 iT F <?3 jT + fT + #5 [T 


(37) 


u 


where R v is the universal gas constant. Since 


_ r 

Jo C " 


dT 


Jj 
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*--r 

U Jn 


'"Pi 


dT 


where s , is the entropy. 
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A. 

Rv 


T - a u + 


n i 


T + 


a 3 i t 2 
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a 4i t 3 a S i ji 

4 5 


a 6i 

T 


(38) 



Q] ,■ In T + a-i iT + — ^ — T ~1" 



+ 



+ 


a 7/ 


(39) 


where ch*Ru, a 7 ,*i?y are, respectively, the enthalpy of formation and entropy at 0 K. 
Other thermodynamic properties can be easily calculated from c Pl , h„ and s,. 

The equations for the transport properties are given in Appendix B. A 
polynomial fit of the logarithm of the property to the logarithm of the temperature is 
chosen for the calculation. They are 

4 

In H, = XX <( ln ^ 

k-\ 

The k, and D tl can be expressed in a similar way. However, since the thermal diffusion 
ratios, 0, are only slightly dependent on temperature, we use polynomials in 
temperature, rather than the logarithm of temperature, i.e. 

4 

In ©i = (41) 

Jk=l 


E. THE COMBUSTION MODEL 

1. One Step Overal Combustion Model . A one-step, forward overall gas-phase 
chemical reaction is assumed for the combustion model. This simple chemical reaction 
system involves a reaction between two reactants (fuel and oxidant) in which they 
combine, in fixed proportions by mass, to produce a unique product. 
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Fuel + Oxidant -» Product (42) 

This model accords with reality with respect to the gross behavior, but suppresses 
the distracting intermediate details. The purpose of so doing is to generate quantitative 

predictions of combustion phenomena which are easy to perform and understand, 

which fit reality in its main features, and which can be refined when necessary. The 
following chemical reaction is considered for the hydrogen diffusion flame under this 
model. 

H 2 + y0 2 -* H 2 0 (43) 

From this model, the production rate w'" in Equation (31) for hydrogen can be 
expressed as 

*H - ~Z H p 2 Y H Y 02 cxp(-E H IR v T ) (44) 

where Z H is the reaction rate constant and E H is the activation energy. Both Z H and 
E h are determined by comparing with experimental data. 

In order to find the source terms for 0 2 and H 2 0 , we introduce the equation of 
conservation of the production rate of all species 


ns 



i—\ 


where y im is the number of atoms of the mth element in the /th species, M, is the 
molecular weight of ith species, / is numbers of element in the reaction, and ns is the 
number of species. 
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From the one step combustion model, we have 


w, 


H, 


W i 


H,0 


M, 


H, 


A/, 


= 0 


H,0 


(46) 


Mh 2 o 


w, 


+ 2 - 


o. 


M 




= o 


(47) 


These will be the production rates for H 2 0 and 0 2 . 

2. Complete Reaction Mechanisms . An eleven reactions model is employed for 
detailed calculations. These reactions have been used successfully in the well-stirred 
reactor combustion simulation [24]. The general form for reversible (or irreversible) 
reactions can be represented by 

ns ns 

2 (m = 1,2,... , nr) (48) 

i=\ f=i 

where o im are stoichiometric coefficients, x> is the chemical symbol for /th species, and 
nr is the total number of reactions. 

The production rate w’" can be written as a summation of the rate of progress 
variables for all nr reactions involving the /th species: 

nr 

*/" = (49) 
m= l 


where u, m = (o",„ - v',„) 

The rate of progress variable, q m for the mth reaction is given by the difference of the 
forward rates and the reverse rates as 
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n£ ns 

q.-^.nw'"-‘..nw"" < 5 °) 

;=i i = i 

where [X,] is the mole concentration of the /th species and k /m and k rm are the forward 
and reverse rate constants of the mth reaction. They are equal to 

k /m = A m T fim exp (-E JR V T) (51) 

and 



where the pre-exponential factor A m , the temperature exponent /?„, and the activation 
energy E„ come from experiment and are given in Table I, P alm denotes atmospheric 
pressure, and k c „ is the equilibrium constant. 

When a third body is needed for the reaction, the rate of progress variable is given 
by the equation below. 



If all species contribute equally as third bodies, then all a im = 1. However, it is often 
the case that some species act more efficiently as third bodies than do others. The <x im 
coefficients are then used to specify the increased efficiency of the /th species in the 


mth reaction. 
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Table I. CHEMICAL REACTION MECHANISMS FOR HYDROGEN-AIR 

DIFFUSION FLAME 



Chemical Reaction 

A m 

0. 

E m 

1. 

H + 0 2 ^ 0 + OH 

5.1E+16 

-0.82 

16510. 

2. 

H 2 + O ^ H + OH 

1.8E+10 

1.0 

8830. 

3. 

H 2 + OH £ H + H 2 0 

1.2E + 09 

1.3 

3630. 

4. 

OH + OH 0 + H 2 0 

6.0E + 08 

1.3 

0. 

5. 

H + OH + M £ H 2 0 + M 

7.5E + 23 

-2.6 

0. 

6. 

H 2 +M £ H + H + M 

2.2E+12 

0.5 

92600. 

7. 

H + 0 2 + M £ HO, + M 

2.1E+ 18 

-1.0 

0. 

8. 

ho 2 + H £ H 2 + 0 2 

2.5E+13 

-0.0 

700. 

9. 

H0 2 + H £ OH + OH 

2.5E+14 

-0.0 

1900. 

10. 

H0 2 + O £ OH + 0 2 

4.8E+13 

-0.0 

1000. 

11. 

ho 2 + oh ^ h 2 o + o 2 

5.0E + 13 

-0.0 

1000. 


For the numerical method implemented in this analysis, the production rates are 
divided into a creation rate and a destruction rate, i.e., 


W; 


p rrr r\ rrr 


(55) 


where 


rrr 


1 j=\ m—\ j— i 


(56) 
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nr nj nr n* 

D"' = 

m = 1 j=\ m=l 7=1 


(57) 


For third body reactions, similar to (54), each sum in the above equation is multiplied 
by the factor 

ns 

f=l 


These production rate calculations will be provided by the CHEMKIN routines. 


3. The Damkohler Number . In the diffusion flame, the chemical reaction is 
generally faster than the species diffusion velocity, that is, the chemical reaction time 
is smaller than the diffusion time. Consequently, the chemical reaction occurs in a 
narrow region between the fuel and the oxidant; the concentrations of the reactants 
are very low in the reaction zone, and the combustion rate is controlled by the rate at 
which fuel and oxidant flow into the reaction zone. 


Here, we introduce the Damkohler number which plays an important role in the 
prediction of extinction condition and in the interpretation of the numerical results. 
The definition for Damkohler number as proposed by Fendell [5] is 


the characteristic flow time 
the characteristic chemical time 


ZnPe 
a(rt + 1) 


(58) 


When D, approaches zero, the flow becomes a frozen flow; if the Damkohler number 
approaches infinity, the flow is in chemical equilibrium. Within this semi-infinite range, 
0 to oo, the actual Damkohler number is determined by the finite chemical rate 
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calculation. Therefore, if the extinction of the diffusion flame occurs at a certain critical 
value of D, for various combinations of fuel and oxidant, the characteristic chemical 
time, and accordingly the chemical rate can be estimated by measuring the flow time 
at flame extinction. 

F. NUMERICAL ANALYSIS 

1. The Collocation Method . The collocation package COLSYS was developed 
for solving mixed-order systems of multipoint boundary- value problems [16]. This is in 
contrast to several other codes which required conversion of a given problem to a first 
order system, thereby increasing the number of equations and changing the algebraic 
structure of the discretized problem. Numerous numerical experiments have 
demonstrated the stability of the collocation procedure, including adaptive mesh 
selection and error estimation. Therefore, this efficient package is chosen to solve the 
momentum equation with its boundary conditions (Eq.33 and 34). 

The method of spline collocation at mesh points, as described in detail in Ref. 17, 
is implemented in COLSYS to solve the nonlinear differential equations. The problem 
is solved on a sequence of meshes with the solutions expressed in terms of a polynomial 
basis. After each iteration, the error is estimated to check against a user-prescribed 
tolerance. If deemed worthwhile, a redistribution of the mesh points is performed to 
roughly equidistribute the error to each subinterval and the solution is recomputed. If 
not, each subinterval is halved, a new solution is computed and the error is 
re-estimated. This iteration procedures will continue for linear problems until the error 
is within the tolerance. For nonlinear problems, the damped Newton's method of 
quasilinearization is performed. Thus, at each iteration a linearized equation is solved 
by the same collocation procedure as described above. 
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Since neither C nor p,/p in Equation (29) varies significantly with mixture 
composition, Equation (29) is partially 

uncoupled from the energy and species equations and it will be solved only after 
certain number of iterations to account for the effect of the similarity function/on the 
solution. 


2. Time- Marching Technique . In order to employ the time dependent technique 
for the coupled set of energy and species equations, Equations (30) and (31) are 
replaced by the time-dependent form 


1 dg 

a dt 


C 


Sr 


+ fSr, + 



{L e -\)Y iv + 


L t Y- 

u e, 1 1 



(59) 


lIXl 

a dt 




+ fY„ + 


\Pr e dv )„ + 


(n + l)a P c i, 


(60) 


Because the concentration and temperature gradients are not uniform in the flow 
field, a continuously collapsing grid is required. The grid spacing is chosen to be fine 
toward the stagnation plane, where the gradients are much steeper. The collapsing 
factor, which is defined by (A*/,*, — A^)/A?/ ; , should never exceed 10 percent to prevent 
the unsymmetrical truncation error accumulating rapidly. The finite difference 
algorithm is based on the control volume approach [25]. The control volume and the 
grid spacing distribution is shown in Fig. 3. 

Equations (59) and (60) can be converted into the finite difference form through 
discretization. We use central differencing and upwind differencing methods for the 
diffusion term and the convection term, respectively, while the system of equations will 
be treated in a fully implicit manner with the exception of the source term for I I 2 0 
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Figure 3. The Grid Spacing Distribution along rj direction, 
when the one-step combustion model is used. The equations can be expressed in the 
following form 


E,- X g£\ + EjgJ +l + E J+]g ;:; = a 



E ’+' “ "“T/^oo] 

Ej - + fi) - E i - 1 - £ y*i + (/» _ A) 


where max[X , Y] denote the greater of X and Y 
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Sij-i = ~ max[ -(S c ) w ,0.0] 

S iJ+l = - maxQ( S c ) e , 0.0 ] 

' ' e 

^ i,J = A ty-l(‘^T fj _ ^l,J— 1 — ^f,y+l C(^c)w — (^c)jll 


where 


S c =/ + 


I T 

C L *. d6 

p r e d v 


D c=fj 

d c=f/ + fy"iiY, 


for H 2 0 in the one-step model 
otherwise 


B i = A 


r'. 

a At 



where 

{ C c = 0.0 for H 2 and 0 2 in the one-step model 
C c = C"’i otherwise 
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Equations (61) and (62) constitute a ( j-2 ) x ( j-2 ) block tridiagonal matrix, 
where each block is a ns x ns matrix. This system will be solved by the Thomas 
Algorithm. 

In order to compare the calculated normal velocity distribution with the 
experimental data, we need to find the physical normal coordinate that corresponds to 
i/. From Equation (19) we have 

dr] _ f p e a(w +1) ) 2 p 
dy ( j Pe 


O' 

f n £*_ A f PA^± 0 1 T P , 

J 0 P ^ j J J 0 ^ (63) 

This gives us the relationship between rj and y. The Romberg method will be used to 
perform the numerical integration of this equation. 

For the global time-marching procedure, we start with initial guess profiles for 
both the temperature and the species mass fractions. A linear distribution with an 
assumed flame front as the peak point for both the profiles is chosen for this purpose. 
These initial profiles will then be smoothed by a commonly used curve-fitting routine. 
In order to find the coefficients appearing in the governing equations, The CIIEMKIN 
routines are implemented and the same curve-fitting routine is used to find their 
derivatives. The momentum equation is, then solved by COLSYS to find the similarity 
parameter f The equation of motion is partially decoupled from the energy and species 
equations, it will be solved after certain numbers of time steps to update the value of 
f The time marching technique is implemented to find the enthalpy and species 
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distributions, which will be followed by a check for steady-state solutions. The 
criterion for testing for steady state solution is as follows. 

| Y‘ +s — Y‘ | < tolerance (64) 

If these values are not within the specified tolerance a new temperature distribution 
will be calculated from the enthalpy using the interval-halving method and the 
coefficients will be updated and the procedure will be continued. Once Equation (64) 
is satisfied, the solution is assumed to have reached steady state and the normal 
velocity v is calculated from Equation (32). The solution obtained for one case will be 
used as the initial guesses for other cases; a series of solutions may be obtained in this 
manner, each time varying the parameters only slightly. This procedure greatly reduces 
computer time and effort involved in running a case. 
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IV. RESULTS AND DISCUSSION 

A. THE EXPERIMENT 


Several Opposed-Jet Burner (OJB) experiments were conducted by Pellett et al. 
[18, 19, 20, 21, 22] at the NASA Langley Research Center. The main experimental 
equipment included two, equal diameter, coaxial, circular tubes. Three different sets 
of OJB tubes were used, 2.7 mm, 5.0 mm and 7.0 mm i.d., respectively, in order to cover 
a wide range of conditions. The two OJB tubes were set in an aluminum mounting 
such that only axial adjustments in tube separation distance were permitted, while the 
radial movements were mechanically restricted. The burner was bathed in an argon 
atmosphere to prevent extraneous combustion outside the central impingement zone, 
which would hinder visibility of the flame. All measurement were made at atmospheric 
pressure and the reactants and tubes were maintained at room temperature throughout 
the measurements to preserve consistency in experimental conditions. A digital mass 
flow metering system was implemented to produce accurate flow rate for both tubes. 
The flow rates of fuel and air were adjusted such that the dish-shaped flame was 
centered between the two tubes. The flow was always laminar. Since the flow through 
the circular tubes is laminar at low Reynolds number, the velocity profile at the tube 
exit is parabolic. 

When a laminar jet of nitrogen-diluted hydrogen was ejected from one tube and 
pure or contaminated air from the opposed tube, a counterflow flame was centered 
between the two tubes. For every specific hydrogen mole fraction measurement, the 
mass flowmeters were used to control this fixed rate. The respective jet flow velocities 
were increased slowly by adjusting the nitrogen-diluted hydrogen and air flows until 
extinction occurred, at which time the volume flow rates of hydrogen, nitrogen, and 
air were recorded. The flow velocities were then reduced in a similar way so that the 
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torus-shaped flame restored back to the jet axis. This abrupt restoration of the central 
flame is the so called RESTORE condition. The flow rate of each component was also 
recorded after RESTORE was established. In order to employ LDV technique to 
measure velocity distributions along the central axis, Al 3 0 2 -seeded air was used. After 
the specific flow condition was stabilized, the normal and radial velocities were 
recorded. 

B. RESULTS 

The present study was undertaken in order to understand the various aspects of 
hydrogen-air combustion such as the flame extinction characteristics. The presence 
of contaminants such as water vapor and carbon dioxide, the role of thermal diffusion 
of light species, and sensitivity to reaction mechanisms are some of the aspects 
considered to be important in this study. The development of an algorithm to solve a 
set of very stiff ordinary differential equations arising from the consideration of detailed 
chemical kinetics associated with the combustion of hydrogen was considered 
important because of the potential for the numerous studies that can be conducted to 
investigate the fundamental mechanisms concerning laminar and turbulent flames. 
These objectives have been accomplished in a two-step process; the first step involved 
the development of a one-step reaction model and the second step involved the 
development of a model with detailed chemistry consisting of 8 species and 22 reactions 
(counting both the forward and the backward reactions). The one-step reaction model 
proved to be very beneficial because it was possible to use the results from it as initial 
profiles for the detailed-chemistry model. Comparisons are made with the results 
obtained from these two models. The results are presented as temperature and species 
profiles and axial velocity profiles. The relative location of the flame with respect to the 
stagnation plane, flame thickness, the maximum values of temperature and species 
mass fractions and their relative locations are some of the items discussed in this 
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chapter. The results presented are only a representative sample to demonstrate the 
capabilities of the computer program rather than a set of detailed case studies. This 
approach was necessary because of the limitations of the size of, and turnaround time 
for the computers available at the University of Missouri-Rolla. The extinction 
conditions were obtained by making several runs and plotting the results on a 
temperature vs. strain rate (or the corresponding Damkohler number) plot. Extinction 
is defined as the condition where the temperature drops precipitously. All the profile 
plots are made with the stagnation point as the origin as shown in Fig 2. The plots 
are made in the physical coordinate rather than the transformed coordinate in order to 
facilitate comparison with experimental data. 

Results from the one-step reaction model and the 8-species model are compared 
in Figs. 4, 5 and 6. Both these calculations were made at conditions corresponding 
flame extinction. Therefore, strain rates of 4500 s~ l and 7700 s ', respectively, were 
used for these two cases. Extinction is seen to take place at nearly the same peak 
temperature for the one-step model and the detailed chemistry model and flame 
thickness is not affected considerably by the choice of the chemistry 7 model. However, 
the strain rate at extinction is higher (7700 .r 1 ) for the 8-species, 22-reaction model 
compared to the extinction strain rate of (4500 s' 1 ) for the single-step reaction model. 
The corresponding extinction strain rate from experiments is 4000 It should be 
emphasized that the activation energy for the one-step model was chosen so that the 

* 

strain rates for the experiment and the model match for one value of the mole fraction. 
Thus, the procedure for the one-step model should be thought of as a means of 
estimating the activation energy for the one-step reaction. The overprediction of the 
extinction strain rate by a factor of two by the 8-spccies model is probably due to the 
fact that the 'effective strain rate' at the flame is higher than the strain rate calculated 
from the cold flow velocity at the tube exit; this point has been discussed by 
Dixon- Lewis [13]. 
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Flame location is an important parameter necessary to discuss flame 
characteristics. We may define flame location as the position of the peak temperature. 
Spalding [3] gave an analytical expression for the flame location from his flame sheet 
theory. Hahn, Wendt and Tyson [12] found it necessary to compromise flame location 
in order to compare experimental results and analysis. In the present study, 
the flame is located approximately 0.015 cm to the right of the stagnation plane, as 
can be seen in Fig. 4. The mass fraction profiles of the major species are shown in Fig. 
5. The one-step model gives a peak H 2 0 mass fraction of about 0.14, whereas, the 
8-species model gives a peak H 2 0 mass fraction of about 0.18. It may be noticed that 
the mass fraction peaks are to the left of the temperature peaks. The mass fraction 
profiles of minor species are shown in Fig. 6. The H-atom profile is seen to be rather 
flat extending from y = 0 to y = 0.03 cm. The H0 2 profile is distinct from the other 
profiles in that it shows two peaks to either side of the flame with their magnitudes 
much smaller compared to other minor species. Hahn, Wendt and Tyson [12] obtained 
analytical results which showed the flame to be thicker than the experimentally 
observed value for methane and CO counterflow diffusion flames. They also observed 
that the experimental reaction zone (as defined by the location of the peak temperature 
and that of intersection of reactant profiles) seemed to be shifted 0.25 cm to the 
fuel-lean side. Several probable factors that might have caused this difference in the 
flame location were identified by these authors; they subsequently adjusted their results 
to make the flame location from experiments and theory coincide. 

The above comparisons indicate that the one-step reaction model gives useful 
information regarding several aspects of the flame. Computer runs with the one-step 
model take much less CPU time compared with the full-chemistry model. The 
full-chemistry model requires very careful preparation of initial profiles for obtaining 
successful steady-state solutions. For these reasons, the results presented hereafter are 
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based on the one-step model. These results could be used as guides for choosing cases 
to be run with the full chemistry model. 

The variation of the maximum temperature as a function of the Damkohler 
number for various values of fuel mole fraction is plotted in Fig. 7. Each data point 
in this figure corresponds to one computer run. The extinction conditions from the 
experiments and analysis are presented in Fig. 8. This plot is obtained from the curves 
shown in Fig. 7. The plot shows the strain rates corresponding to various hydrogen 
mole fractions in the fuel stream. As stated earlier, the activation energy used in these 
calculations using the one-step reaction was chosen so that the results from 
experiments and analysis agree at hydrogen mole fraction equal to 0.6. This procedure 
was necessary because the activation energy for the one-step reaction model could not 
be determined satisfactorily using theory. It is worth mentioning that Tien, Singhal, 
Harrold and Prahl [11] have used a similar procedure to calculate extinction in the 
stagnation point boundary-layer of a condensed fuel using a one-step reaction model. 
The analytical curve in Fig. 7 has a larger slope compared to the experimental curve. 
The underlying reason for this behavior may be traced to the fact that the strain rates 
in this plot are obtained from cold flow at the pipe exit; the 'effective' strain rates may 
be larger than the ones obtained in this manner. The effective strain rate is likely to 
be a function of the heat release and, therefore, also the mole fraction of the fuel. This 
may be the reason why the analysis overpredicts the strain rates at higher mole 
fractions and underpredicts them at smaller mole fractions. It is, therefore, possible 
that improvements in the predictions could be obtained if the strain rates are corrected 
to account for the changes in the velocity profile due to heat release. 





39 



Figure 8. Hxtinction Limits Comparison between I'xpcrimcntul Data and One-Step Model. 
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Figs. 9, 10 and 1 1 show temperature profiles and mass fraction profiles for two 
extreme cases, one near equilibrium and the other near extinction. The equilibrium 
case has a thicker reaction zone because of the increased heat release. It should also 
be noted that the maximum temperature in this case is much higher that the near 
extinction case. Another feature to be observed in these figures is that the flame moves 
towards the stagnation plane as it approaches extinction. This sensitivity of the flame 
location with respect to the stagnation plane to the conditions under which the flame 
exists has been observed also by Dixon-Lewis et al. [13]. They observed that, in this 
regard, agreement for methane-air flame is much better than for hydrogen air flame. 
Pdlett et al. [18-22] observed in their experiments that the flame becomes more planar 
and moves towards the fuel side as extinction is approached. This phenomenon is still 
not properly understood and no simple explanation exists for this behavior. It appears 
that the higher heat release in the hydrogen-air flame as compared with the 
methane-air flame does play a role in making the effect more pronounced in the former. 
Simply shifting the curve in order to match with experimental data, as was done in Ref. 
12, does not appear to be satisfactory. Hahn et al.[12] attributed the discrepancy in 
flame location to factors such as inaccurate calculation of transport properties and/or 
kinetic parameters, incomplete chemical reaction mechanism, and the absence of 
buoyancy and radiation terms in the analytical model. The present results along with 
the observations of Dixon-Lewis et al. and Pcllett et al. indicate that this phenomenon 
requires more careful attention for it to be completely explained. 
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One concern that has arisen in the past regarding the use of hydrogen as the fuel 
is the role that thermal diffusion of light species such as H 2 and H has on the flame 
characteristics. A case study was, therefore, conducted to determine the extent to 
which thermal diffusion affects flame characteristics. The results are plotted in Figs. 
12 and 13. These two figures demonstrate that thermal diffusion has little effect on the 
flame; its only influence appears to be a slight reduction in the peak temperature. It 
appears that using a more detailed model for thermal diffusion than the one presently 
used will not have any noticeable effects on the flame. 

The experimental temperature profiles of Pellett et al. are compared with the 
predictions in Figs. 14 and 15. The experiments are based CARS (coherent antistokes 
Raman sepctroscopy) measurements. These figures and other results not shown here 
indicate the tendency to overpredict the temperature compared to experimental values. 
Hahn, Wendt and Tyson [12] and Dixon- Lewis et al. [13] also have reported 
overprediction of the temperature. Both these groups attributed the differences to the 
fact that no radiation correction was applied to the temperature measurements using 
thermocouples. However, the present results indicate that the discrepancies between 
experiments and analysis cannot be fully explained as being caused by thermal 
radiation because the temperature data used in the present study are based on CARS 
measurements. More work on the influence of grid distribution, number of species 
involved and the reaction set used for the calculations must be performed to answer 
this question. No direct measurement of the flame location with respect to the 
stagnation point was made in these experiments. Therefore, the differences in the flame 
location seen in Figs. 14 and 15 may not be realistic; the velocity profile plots give a 
more accurate picture of the flame location. 
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Finally, the experimental velocity profiles of Pellett et al. [22] are compared with 
the predictions in Figs. 16 and 17. The experimental data are based on LDV 
(laser-doppler velocimetry) measurements. The main feature to be observed in these 
two figures is the presence of the velocity maxima caused by heat release in the flame 
zone. Differences in the flame location between the experiments and the calculations 
persist in these figures. There is fairly good agreement in the magnitudes of the peak 
velocity. 
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V. CONCLUSIONS 

A capability has been developed to analyze counterflow diffusion flames. The 
ability to solve the stiff equations governing this problem when hydrogen is used as the 
fuel, is considered to be an important step in understanding the complex phenomena 
associated with counterflow diffusion flames. The numerical procedure based on the 
control volume approach with time-marching appears to be better suited to the present 
problem than the methods previously used for the solution of stiff ODEs. The present 
model uses detailed chemistry and accounts for variation of Prandtl number and Lewis 
number as well as considers the effect of thermal diffusion on the flame. It is shown 
that a one-step model can predict several features of the flame, while the detailed 
chemistry model can be used for fine tuning the results. The present results show that 
thermal diffusion has negligible effect on the characteristics of the flame.. Further 
investigation is needed to settle some issues such as the flame location with respect to 
the stagnation plane, overprediction of the temperature by the analytical model and the 
considerable disagreement between theory and experiments on the extinction 
conditions; it is suggested that the actual strain rate in the flame zone, which is higher 
than the cold flow strain rate, must be the one to be considered for analyzing 
extinction. 

Several valuable studies may be conducted using the present analytical model. 
The effect of contaminants on the flame is one such study currently being done. 
Consideration of more species such as H2O2 and more reactions could be used to 
understand the effect of the chemistry model on the flame. Use of the present model 
to investigate a recently proposed model for turbulent diffusion flames, according to 
which turbulent diffusion flames consist of a cluster of strained, laminar diffusion 
flamelets, will be an interesting and challenging task for investigators interested in 


turbulent diffusion flames. 
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APPENDIX A 


DERIVATION OF THE SIMILARITY TRANSFORMATION EQUATIONS 
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We will outline here the detailed derivation of Equation (29) through (32). For 
convenience, we rewrite the operators (21), (22), and introduce a more general 
expression for stream function 


_d_ = put*" d 
dy (2£) 1/2 d V 


( 21 ) 


_d_ 

dx 


2n 

PePiPe* 


d d 

d{ dr\ dt, 


( 22 ) 


* = (2£)' /2 /(£, V) 


(A. 1) 


A. NORMAL VELOCITY 

The equation for normal velocity is 

* d* 
pvx = - — 

From Equation (22) and uju e = d/jdrj, we have 

= s/?4 + fr, Vi) + PePeUe* 2 " f!s[U 

= P e Pe u e x2n (J2 T/ ( + -fi {"/,». +n-M) 


Then, from Equation (25) 


pv 


-n jj 

X dx 

- PeW \Ju a + -Mf.Vi +//V u) 


(25) 


(A. 2) 


(A. 3) 



is the case, we neglect those derivatives with respect to £ 


PV = - PePeUgX" flyfit 


( 32 ) 


B. MOMENTUM EQUATION 


The equation of motion is 


puu x + pvu y = ~ Px + (pu y ) y 


( 2 ) 


From Equation (21) and (22), we have 

■Jj = PePe U W n + ftp, + PePe^Vr, \ ( /I J ) 


2 n 

du_ = Pjfe± f 

dy Jt " 


(,1.5) 


d(pUy) pUe x n ( u]x n r ^ 

* “ Ju V" M f ”r 

3 2/i r j ,\ 

P“ e * / PM r \ (/ 16 ) 

3 2 * 

= PeMe 2( e ( C /^)r, 


The first term on LHS of Equation (2) becomes 

puu x = pu e f n [p e p e u 2 x 2n (/ ?n + f vrj rj { ) + p e p e u e x 2n f r] « e J 


(/1.7) 
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Combining Equation (A. 3) and (A. 5), the second term on the LHS is 

pVUy = - p e p e pu e x 2n (f i +//2£ + fr,Vt )fr,r, (^- 8 ) 


dp_ 

dx 


du e 

- p ‘ u '-d7 


2 2 In 

= ~ PeHe u e x 


<K 

dt, 


(A.9) 


Dividing each term by p.p.pulx^jl^, we get the following momentum equation 

(C/„ (ir - /, 2 ) “ - /;/,,) 1°) 

For / =f{y) and u, = ax, we have the final momentum equation. 

(C/„), + //„ + - /, 2 ) = 0 (29) 


C. ENERGY EQUATION 


The energy equation (14) is 


+ 4 )* + p, ( h + -rl 

= {(Klc~ p ){h + u 2 /2),} + {(l/2)[> - (*/<p]( U 2 )^ 

+ Ap - (*/cp](A + 

+ 


(14) 


61 


We can have following relations from Equation (21) and (22) 


djh + m 2 /2) 
dx 


= Pe^e u eX 2n h e ( g { + 


(-4-H) 


d{h + u 2 I2) pugx n t 

*y Ju eS " 


(A. 12) 


8c i = p^y_ y 

dy Ju Ci ' l '' 


(/U3) 


dT 
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PU<r x 

■M 


T e d„ 
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(,1.14) 


S(u) 

d y 


8u 2 r P U e^ r 

= 2 “77 = 2u - f --^r f ’ 




(,1.15) 


The first term on LHS becomes 


pu{h + u j 2)^ p e p e pu g x T &rj %) 


(,1.16) 


With Equation (A. 3), the second term on LHS is 


H* + *0'2), 

- - p tlie uy(j2( f, + + flj2i)^-h,g, 

PePeP^e^ ~E S^f "E frfg^^l^) 
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Similarly, the terms on the RHS become 
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Dividing each term by p,p t pulx in h.l2£, we get the following energy equation 

( 'a*’)* +/s ’ + - 1)! % + ~ir' f’jj, (/)22) 

+ (»X){( i - pr’)c/x}, - «U* { -/[*,) 


If we assume wj « /i, and neglect derivatives with respect to we have Equation 


(30). 



+ /&, + 



£^1^0 

0 


= 0 (30) 

>1 


D. SPECIES EQUATIONS 


The boundary layer species equation is 

puc ix + pvCjy = | Dfic iy + D[pc t T y jT^y + w{” 


(15) 
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From Equation (21), we find 


= PeHe“eX 2nc i,e{Yii + Y it , >?{) 


The first term on LHS of Equation (15) becomes 

P UC ix = P« e /,P e w\ e (^ + Y ir,Vi) 

= PePeP u W n Ci,e( Y iS + Y ir , > 7 |) 

Combining Equation (A. 13) and (A. 5), the second term on the LHS is 
pvCjy p e fl e pUgX ~b //2£ T yjj 

The first term on the RHS becomes 
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04.26) 


Dividing each term by p e p e pnJjc 2n c,,,/2^ , we get 
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jC 

Pr 



+ /Y, + 


( c de \ , 

[p r e d V Jr, + 


(rt+l)a P c i,e 


= 2{(/,4 - /;/„) 


If we neglect derivatives with respect to , we have Equation (31). 
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This complete all the derivations. 
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THE TRANSPORT EQUATIONS 
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A. VISCOSITY 


The viscosity of species / is given by the expression [26] 


Mi = 


5 

16 


(5.1) 


where o, is the Lennard-Jones collision diameter, and W, is the mass of the molecule. 
The collision integral depends on the reduced temperature given by 



KbI 

£ / 


(5.2) 


ard the reduced dipole moment given by 


where r, is the dipole moment, £, is the Lennard-Jones potential well depth, and K B is 
the Boltzmann constant. The integral value will be given by the quadradic interpolation 
of the previous experimental work by Monchick and Mason [27], 


B. THERMAL CONDUCTIVITY 


The thermal conductivities then follow from the expression [28] 


K i ^ (^trans c v,lrans "b 5 ro[ C vrol + Z\.,^ C vv jt)) 


(5.4) 


where 
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(5.5) 
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1 rot 
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(5.6) 
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Vi 


(5.7) 


and 


r = 


5 
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p£n 

Vi 


(5.8) 


T = U, 


rot 



L v,rot 


3 R 


u 



(5.9) 


The constant volume specific heat c, relationships due to different excitation 
energy are dependent on whether the molecule is linear or not. For linear molecule, 
we have 


c v,rot _ _2 

c v,trans 3 


(5.10) 


^v.rot 

R 


(5.11) 


^vib 



(5.12) 


where c „ is the full constant volume specific heat of the molecule. In the case of a 
nonlinear molecule, we have 


c v,rot 

c v,trans 


(5.13) 
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c v,rol 3 

"XT = 7 


(5.14) 


c v,vib c v 


(5.15) 


For single atoms (like H,0 atoms), there are no molecule energy contribution due to 
rotation and vibration. Hence 


Ki = 



(5.16) 


where F xnni = 5/2. For selF-difFusion coefficient, we evaluate the following expression 
[26] 


3 ^nK^IW, 
16 p no]d x ' *)* 


(5.17) 


The rotational relaxation collision number t/ ro , is dependent on temperature. If 
data is available at 298 K, we have 


U tol (T) = 


^ro t (298) 


A(298) 

"A(7T 


(5.18) 


where 
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+ ri 


3/2 


f/K £ 

T 
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C. DIFFUSION COEFFICIENTS 


The binary diffusion coefficients can be evaluated from the following expression 


[26] 


3 

16 2 rvO.i)* 

pna ki Q K 


(5. 20) 


where W kl is the reduced molecular mass for the (k, i) species pair 

2 W k W , 

Wu = — — 

kl W k + 


(B- 21 ) 


and a kl is the averaged collision diameter. In computing the averaged quantities, we 
consider two cases, depending on whether the collision partners are polar or nonpolar. 
For the case that partners are either both polar or both nonpolar, we assume 


£ k£ 

K s 



(5.22) 


a ki ~ 
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(5.23) 
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(5.24) 


For the case of a polar molecule interacting with a nonpolar molecule, we have 
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np 


K B 



(5.25) 
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(5.27) 


where 


c 


1 + 


7“» t p 



(5.28) 


In the above expression is the reduced polarizability for the nonpolar molecule, and 
Tp is the reduced dipole moment for the polar molecule. The reduced values are given 
by 


(5.29) 


V £ A 


(5.30) 


The following reduced parameters are required for the collision integral Q 0 - 11 ’ 
evaluation 



£ ki 


(5.31) 


S* _ 1 * 2 

°ki ~ 2 


(5.32) 


Now we introduce a correction factor to the binary diffusion coefficients [29|. 
Equation (B.4) becomes 

D ki = D ki (1 + X kl ) (5.33) 

where the correction factor is given by 
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D. THERMAL DIFFUSION RATIOS 


The thermal diffusion ratios can be defined by introducing the thermal diffusion 
velocity 


Di®i 1 3T 
X t T dx y 


( 5 . 42 ) 


where x v is a spatial coordinate and X, is the mole fractions. We only consider the 
thermal diffusion for the species / whose molecular weight is less than 5. The thermal 


diffusion ratio is given by [30] 


(S.43) 

k±i 


where 


fik i 


15 + 5)(6Q f — 5) M k — 

2 I2^ 1 + 55) 


(8.44) 


These are the semi-experimental approximations used in the transport properties 


evaluation. 



